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ABSTRACT 


A review  of  Mathematical  Function  Library  for  Microsof t-FORTRAN , John 
Wiley  & Sons,  1989,  and  associated  computer  software  is  presented.  The  pack- 
age consists  of  xvii  + 341  pages,  25  1/2  cm,  loose  3-hole  punched  leaflets  in 
a ring  binder  and  three  5 1/4"  diskettes.  Price:  $295. 
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DISCLAIMER 

This  paper  presents  the  results  of  research  which  produces  a 
preliminary  evaluation  of  numerical  computing  software  made 
available  in  a commercial  product.  The  research  is  published  to 
coordinate  and  encourage  work  in  the  important  area  of  software 
quality;  the  results  of  the  research  are  not  to  be  interpreted  as 
an  official  assessment  by  the  NIST. 
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1. 


INTRODUCTION  AND  DESCRIPTION  OF  THE  UL  LIBRARY 


In  1964  the  National  Bureau  of  Standards  (now  the  National  Institute  of 
Standards  and  Technology)  issued  a massive  handbook  of  formulas,  graphs  and 
numerical  tables  of  the  elementary  mathematical  functions  and  the  so-called 
higher  transcendental  functions  or  special  functions  of  mathematical  physics 
[6].  This  NBS  Handbook  immediately  filled,  and  continues  to  fill,  a tremen- 
dous need  in  scientific  work;  according  to  the  Science  Citation  Index  (pub- 
lished by  the  Institute  for  Scientific  Information,  Inc.,  Philadelphia,  PA), 
the  current  rate  at  which  it  is  cited  in  the  mathematical  and  scientific  lit- 
erature is  of  the  order  of  1,300  entries  per  year. 

The  need  for  the  numerical  tables  of  the  elementary  functions  in  the  NBS 
Handbook  has  now  largely  disappeared:  library  routines  for  generating  expo- 
nentials, logarithms,  and  trigonometric  and  hyperbolic  functions  are  available 
in  all  major  scientific  software  libraries  as  well  as  being  required  by  the 
FORTRAN  Standard  [1].  Also,  these  routines  are  incorporated  in  hand-held  cal- 
culators designed  for  scientific  calculations. 

The  loose-leaf  manual  and  accompanying  diskettes  under  review,  which  we 
shall  refer  to  as  the  "UL  Library",  may  be  regarded  as  an  attempt  to  replace 
the  numerical  tables  of  the  higher  transcendental  functions  supplied  in  the 
NBS  Handbook  by  a comprehensive  software  package.  The  functions  treated 
include  Bessel  and  related  functions,  hypergeometric  and  confluent  hypergeo- 
metric functions,  elliptic  functions  and  integrals,  exponential  integral  and 
related  functions,  error  function  and  related  functions.  Gamma  and  incomplete 
Gamma  functions,  orthogonal  polynomials,  probability  functions  and  random  num- 
ber generators.  This  list  is  not  identical  to  the  list  of  functions  tabulated 
in  the  NBS  Handbook;  among  the  omissions  are  parabolic  cylinder  functions, 
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Mathieu  functions,  spheroidal  wave  functions  and  the  Riemann  Zeta  function. 

The  UL  Library  is  designed  to  be  used  on  personal  computers  of  IBM  type, 
equipped  with  a Microsoft  FORTRAN  77  compiler.  For  efficiency,  a numeric  co- 
processor is  recommended.^  The  operating  precision  is  IEEE  double  precision 
(53  bits  in  the  floating-point  mantissa),  but  the  accuracy  of  the  computed 
function  values  is  generally  less,  sometimes  well  below  single  precision. 

The  UL  library  is  an  extremely  ambitious  project,  especially  as  it  ap- 
pears that  all  of  the  programs  have  been  constructed  from  scratch.  For  such  a 
project  to  be  completed  successfully  its  authors  need  to  have  a thorough 
knowledge  of,  and  experience  in,  several  areas  of  classical  and  numerical 
analysis,  including  analytic  properties  of  the  higher  transcendental  func- 
tions, asymptotic  analysis,  approximation  theory,  and  error  and  stability 
analyses.  How  well  have  the  present  authors  (who  are  nameless)  succeeded? 

2.  NEGATIVE  RESULTS  OF  NUMERICAL  TESTS 

A comprehensive  answer  to  the  question  Just  posed  would  necessitate  a 
tremedous  amount  of  numerical  testing.  We  concentrated  our  testing  on  Just  a 
few  functions  with  which  we  have  had  previous  software  experience,  namely 
Airy,  Bessel,  hypergeometric , confluent  hypergeometric  and  Legendre  functions. 
Usually  the  UL  Library  performed  in  accordance  with  the  specifications  for 
each  routine.  However,  we  studied  the  documentation  in  detail,  looking  for 
major  ways  in  which  the  algorithms  used  might  go  astray.  Unfortunately,  we 
were  successful. 


The  library  diskettes  provide  Fortran  programs,  which  could  be  com- 
piled and  used  on  a personal  computer  without  a co-processor,  as  well  as  co- 
processor assembly  code  for  each  subroutine.  For  the  purposes  of  this  review, 
attention  is  restricted  to  the  co-processor  assembly  code. 
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2. 


THE  FIRST  TYPE  OF  FAILURE 


The  first,  and  very  serious,  type  of  failure  arises  when  the  UL  Library 
delivers  completely  wrong  answers  without  any  error  message.  An  example  oc- 
curs with  the  routine  BESJR  in  §8.8,  which  is  designed  to  generate  the  Bessel 
function  real  values  of  the  argument  x and  order  v.  With  x = 

971  and  v = 2^,  4^,  6^,  8^,  10^  BESJR  yields  the  following  values  of  J^(x), 
to  ten  decimal  places: 

-0.00670  05817,  0.03764  71379,  -0.05418  00781,  0.31114  25305,  -0.06567  89923. 
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These  should  be  compared  with  the  correct  values 

0.01592  10880,  -0.05237  32560,  0.10316  94874,  -0.14727  18330,  0.14363  19977. 

Not  only  are  the  numerical  values  totally  incorrect,  all  the  signs  are  wrong, 
too.  Similar  gross  errors  occur  for  v = 22^(2)50^.  On  the  other  hand,  BESJR 
generates  correct  values  (within  the  prescribed  error  tolerance)  for  the 
intermediate  values  v = 12^(2)20^,  and  also  for  ^ ^ ^ 

again  being  9n. 

The  reason  for  these  errors  appears  to  be  that  J.  C.  P.  Miller’s  back- 
ward recurrence  algorithm  has  been  used,  with  the  trial  values  normalized  on 

the  value  of  J (97r).  Since  J (97i)  is  zero,  this  procedure  is  bound  to 
1/2  1/2 

lead  to  meaningless  answers.  Yet  this  cannot  be  the  entire  explanation, 
otherwise  all  values  in  the  range  v = 21^(1)50^  would  be  incorrect  and  not 
merely  alternate  ones.  Thus  the  documentation  must  be  in  error,  too.  And 
there  may  be  further  inaccuracies  here.  For  example,  according  to  the  docu- 
mentation the  value  for  v = 2^  is  computed  from  the  power-series  expansion 
2 

Obtained  by  use  of  D.  E.  Amos’  package  [2] . See  also  comments  made 
below  on  validation. 
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of  J^(x):  either  the  Miller  algorithm  was  used  instead,  or,  perhaps  less 
likely,  the  power  series  was  summed  incorrectly. 

Similar  errors  occur  for  other  values  of  v,  both  integer  and  noninte- 
ger. For  example,  with  x = 8.65372  79129  ...  (the  third  positive  zero  of 
J^(x}),  BESJR  generates  accurate  values  for  v = 0, 1,2,3  and  7(1)50,  but 
grossly  inaccurate  values  for  i>  = 4,5,6.  Furthermore,  erroneous  values  of 
J^(x)  are  generated  when  the  value  of  x is  merely  moderately  close  to  one 
of  the  critical  values,  the  magnitudes  of  the  errors  being  inversely  propor- 
tional to  the  distances  of  x from  the  critical  value. 

A companion  routine  to  BESJR  is  BESYR  (§8.10),  which  is  designed  to 
generate  the  Bessel  function  Y^(x)  for  real  values  of  x and  v.  Since  one 
of  the  algorithms  used  in  BESYR  draws  upon  values  of  J^(x),  we  expected — and 

found — some  difficulties.  For  example,  BESYR  computed  Y^(9n)  with  wrong 

11  11 

signs  and  incorrect  numerical  values  for  = -502(2)-22^  and  -102(2)-2^, 
while  the  results  were  correct  at  intermediate  values  of  v.  (This  similarity 
of  behavior  with  BESJR  reflects  the  identity  Y (x)  = (-)'^J  (x)  when  = 

^ V -V 

n - 2>  being  an  integer.  ) 

2.2.  The  Second  Type  of  Failure 

A second  type  of  failure  is  the  generation  of  results  which,  while  not 
completely  inaccurate,  contain  errors  greatly  in  excess  of  the  accuracy 
claimed  in  the  documentation.  The  routine  BESYR  illustrates  this  point.  The 
documentation  claims  at  least  12-digit  accuracy  when  the  order  v is  an  inte- 
ger. This  claim  is  careless  because  it  takes  no  account  of  the  inevitable 
loss  of  relative  precision  in  the  neighborhoods  of  zeros.  But  the  situation 
is  actually  much  worse.  With  v = 119  and  120,  we  found  that  in  the  range 
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X = 848(0.1)852  the  accuracy  often  fell  to  8 or  9 digits,  even  well  away  from 

the  zeros  of  Y (x)  and  Y (x). 

119  120 

Another  example  is  provided  by  the  routines  AIRYA  and  AIRYAD  (§§8.1  and 
8.3)  for  the  Airy  function  Ai(x)  and  its  first  derivative.  In  both  routines 
at  least  13-digit  accuracy  is  claimed  when  -10^°  ^ x s 100.  But  we  found 
many  values  of  x in  the  (zero-free)  range  5.2  to  5.8  for  which  they 
yielded  only  8 or  9 correct  digits. 

2.2.  The  Third  Type  of  Failure 

The  third  type  of  failure  arises  when  the  UL  Library  generates  an  innac- 
curate  value  and  the  user  is  warned  by  an  error  message  such  as  "unable  to 
compute  the  . . . function  with  acceptable  accuracy"  or  "numeric  overflow  in 
the  . . . function".  If  these  failures  occur  too  frequently,  then  there  will 
be  huge  gaps  in  the  effective  ranges  claimed  for  the  variables.  Such  is  the 
case  with  the  routines  CHGFU  (§9.2)  for  generating  the  confluent  hypergeome- 
tric function  U(a,b,x),  and  HPRGMT  (§20.1)  for  generating  the  hypergeometric 
function  F(a,b;c;x). 

CHGFU  employs  two  algorithms.  The  first,  evaluation  of  the  asymptotic 
expansion  of  U(a,b,x)  for  large  x,  is  quite  sound.  The  second  is  based  on 
a formula  that  expresses  U(a,b,x)  as  a difference  of  two  M-type  confluent 
hypergeometric  functions,  is  unsound  because  of  the  potential  for  massive 
numerical  cancellation.  As  a result,  although  the  documentation  claims  that 
the  effective  ranges  of  the  variables  are  given  by 

-50  :£  a :s  50,  -50  i b 5 50,  -100  :£  x :£  100, 

with  the  exclusion  of  integer  values  of  b and  nonpositive  integer  values  of 
a,  extensive  regions  are  inadmissible.  These  include: 
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a = 0.5,  b = 0.5,  7.4  :£  X :£  19.9;  a = 2.7,  b = 5.4,  10.1  ^ x :s  437.4; 
a = 20,  b = 2.5,  0.34:£x2  20,000. 


For  HPRGMT  the  documentation  states  that  a and  b can  have  any  real 
value  between  -10^°  and  10^°,  c can  have  any  real  value  between  -10^° 
and  10  (other  than  a negative  integer)  and  x can  have  any  real  value 
between  -10^°  and  1.  Again,  these  claims  are  misleading.  For  example, 
when  X is  in  the  interval  [0,1)  HPRGMT  sums  the  hypergeometric  series 


F{a, b; c; x) 


00 


= I 


a (a+1 ) 


♦ (a+n-1 )b(b+l ) * » ♦ (b+n-1 ) 
c (c+1 ) • • • (c+n-1 ) 


n 

X 

n! 


Since  the  radius  of  convergence  is  1 the  algorithm  must  fail  for  values  of 
X sufficiently  close  to  1.  Sample  intervals  of  failure  were  found  to  be: 


a = b = 100, 

c = 1, 

0.95  :s  X < 1; 

a = b = 10000 

c = 1, 

0.0013  s X < 1; 

. . ...8 

-If) 

a = b = 10  , 

c = 1, 

V 

X 

VI 

1 

o 

Even  when  failure  does  not  occur,  execution  can  be  extremely  slow.  For  a = b 
= c = 1 and  z = 0.99999,  20  minutes  elapsed  on  a 25  mhz  IBM  PS2  Model  80 

before  the  answer  was  produced. 

3.  REMARKS  ON  THE  ALGORITHMS  USED 

The  weaknesses  in  the  algorithms  used  for  J^(x),  Y^(x)  and  U(a,b,x) 

are  avoidable.  Robust  software  for  generating  these  functions  is  already 
available;  see,  for  example,  [3],  [8].  Instead  of  normalizing  the  trial  val- 
ues of  J (x)  obtained  in  the  Miller  algorithm  via  a single  value  of  J (x), 

V ° ° i; 

the  identity 
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00 


I 


(i^+2k)r(i^+k) 

k! 


k=0 


[6],  eq.  (9.1.87),  could  have  been  used.  This  is  the  appropriate  generaliza- 
tion of  the  identity 


1 = J (x)  + 2J  (x)  + 2J  (x)  + 


2 
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0 


that  the  authors  use  in  a routine  BESJ  (§8.7)  for  generating  functions  of 
integer  order.  A stable  way  of  generating  U(a,b,x)  is  backward  integration 
of  the  confluent  hyper geometric  equation,  with  initial  values  derived  from  the 
asymptotic  expansions  of  U(a,b,x)  and  5U(a,b,x)/5x  for  large  x. 

In  addition  to  occasional  poor  choices  of  algorithm,  instances  of  poor 
choices  of  the  actual  functions  being  generated  were  observed.  Thus, 
functions  that  exhibit  exponential  growth  or  decay  when  the  argument  x is 
large  would  have  been  better  replaced  by  their  logarithms.  This  would  greatly 
increase  the  threshold  at  which  overflow  or  underflow  occurs.  Such 
functions  include  the  Gamma  function,  r(x),  the  exponential  integral, 

Ei(x),  the  complementary  error  function,  erfc  x,  the  Airy  functions,  Ai(x) 
and  Bi(x),  the  incomplete  Gamma  function  r(a,x),  the  modified  Bessel 
functions  confluent  hype rgeome trie  function 

M(a,b,x).  In  the  case  of  r(x),  a separate  routine  is  given  for  £n|r(x)| 
in  §13.3,  but  since  this  is  constructed  simply  by  taking  logarithms  of  the 
values  obtained  by  the  library  routine  for  r(x),  there  is  no  increase  in  the 
overflow  threshold.  The  logarithm  should  have  been  generated  first! 

Troublesome  singularities  can  sometimes  be  avoided  by  introduction  of 
appropriate  factors.  Each  of  the  functions  M(a,c,x)  and  F(a,b;c;x)  has 
poles  at  c = 0,-1, -2 but  both  M(a, c,x)/r(c)  and  F(a,b;c;x)/r(c)  are 
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entire  functions  of  c;  this  was  pointed  out  many  years  ago  in  [7].  The 
statement  in  §9.2  that  U(a,b,x)  is  not  well  defined  when  a and  b are 
negative  integers  and  |b|  ^ |a|  is  incorrect.  As  noted  in  [7],  p.  258, 
U(a,b,x)  is  entire  as  a function  of  a and  b.  The  poorly  chosen  algorithm 
used  to  compute  U(a,b,x)  may  fail  when  a or  b is  an  integer  or  close  to  an 
integer . 

4.  REMARKS  ON  VALIDATION 

We  could  continue  in  this  vein  and  we  could  also  provide  a substantial 

list  of  typographical  errors  in  the  manual  and  operational  defects  in  the 
3 

software.  Instead,  let  us  turn  now  to  the  topic  of  validation.  Several 
articles  have  been  written  on  the  difficulty  of  checking  software  written  for 
the  generation  of  mathematical  functions;  see,  for  example,  [3].  How  were  the 
UL  Library  routines  validated  by  the  authors?  The  only  clue  supplied  in  the 
manual  appears  to  be  the  statement  on  p.  36  that  the  least  accurate  results  in 
the  output  of  an  algorithm  always  occur  at  the  interface  with  another  algo- 
rithm. Presumably  this  means  that  there  has  been  substantial  cross-checking 
at  these  interfaces.  This  is  indeed  a powerful  type  of  check,  but  one  that 
does  not  guard  against  all  types  of  algorithmic  and  programming  errors,  as  we 
have  already  observed  with  the  routines  BESJR  and  BESYR. 

There  are  two  subsections  (§§6.2,  6.3)  in  which  the  authors  encourage 
users  to  check  output  by  using  identities  satisfied  by  the  higher  transcen- 
dental functions.  For  example,  the  error  and  Bessel  functions  are  both 
special  cases  of  confluent  hyper geometric  functions.  However,  these  kinds  of 

3 

For  example,  there  is  a statement  on  p.  24  that  the  computer  will  halt 
when  a library  routine  is  called  and  a co-processor  is  not  present.  We  found 
this  is  not  always  true;  furthermore,  instead  of  identifying  the  absence  of  a 
co-processor  as  the  problem,  the  error  messages  misleadingly  report  errors  in 
the  library  routine. 
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identities  are  rather  specialized,  and  there  is  always  the  danger  that  they 
may  have  already  been  used  in  constructing  the  library  routine.  For  example, 
there  are  identities  that  relate  the  Airy  functions  Ai(x),  Bi(x)  and  their 
derivatives  to  Bessel  functions  and  modified  Bessel  functions  of  orders  ±1/3, 
±2/3.  But  if  these  identities  are  used  as  cross-checks  then  the  inaccuracies 
we  noted  above  for  the  routines  AIRYA  and  AIRYAD  may  not  show  up  because  the 
same  kind  of  inaccuracy  is  present  in  the  (parent)  routine  BESKR  (§8.14)  for 
the  modified  Bessel  function  K^(x).  In  contrast,  a powerful  form  of  cross- 
check that  is  not  mentioned  in  the  manual,  but  which  is  widely  applicable,  is 
to  employ  identities  of  Wronskian  or  Casoratian  type.  Examples  are: 

Ai(x)Bi'(x)  - Ai'(x)Bi(x)  = I/tt, 

J (x)Y  (x)  -J  (x)Y  (x)  = 2/(7rx). 

u+i  V u v+i 

Checks  of  this  type  were  used  to  detect  errors  in  the  UL  Library  routines. 

They  were  also  used  to  validate  results  from  other  libraries  used  to  compare 
against  the  UL  Library. 

5.  SUMMARY  AND  RECXIMMENDATIONS 

In  summary,  the  UL  Library  provides  PC  users  with  a new  set  of  routines 
for  generating  an  extensive  collection  of  higher  transcendental  functions, 
indeed  most  of  the  functions  tabulated  in  the  NBS  Handbook.  The  library  is 
relatively  easy  to  install  and  its  price  is  reasonable.  It  will  provide  a 
useful  tool  for  mathematical  physicists  and  other  scientists,  especially  those 
engaged  in  calculations  of  exploratory  type. 

However,  in  a comparatively  small  sample  failures  of  various  kinds  were 
encountered,  including  some  extremely  serious  ones.  It  may  be  that  the 
authors  lacked  some  experience  and  expertise,  or  perhaps  just  the  necessary 
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time,  for  the  mammoth  task  of  constructing  and  testing  a robust  library  of  the 
kind  intended.  Therefore,  users  will  need  to  exercise  great  care  with  any 
output  from  the  library,  applying  independent  checks  wherever  possible. 
Furthermore,  users  must  also  be  prepared  for  disappointments:  the  viable 
ranges  of  a routine  may  turn  out  to  be  a good  deal  less  than  is  claimed  in  the 
documentation,  especially  in  the  case  of  functions  that  have  not  been  treated 
by  earlier  software  workers. 

For  heavy  systematic  computations  many  users  will  find  the  more  robust 
IMSL  and  NAG  libraries  [4],  [5]  to  be  preferable.  The  variety  of  functions 
covered  in  these  libraries  is  not  as  large,  but  the  viable  ranges  of  the 
variables  are  considerably  more  extensive  and  the  precision  is  often  higher. 
Moreover,  both  IMSL  and  NAG  provide  many  desirable  features,  such  as  linear 
algebra  packages,  in  addition  to  routines  for  generating  higher  transcendental 
functions. 
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